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Abstract 

This paper deals with a new fictitious domain coupling algorithm between a 
rigid body and an unsteady compressible fluid flow. The coupling with a rigid 
body is a first step towards the coupling with a Discrete Element method. The 
flow is computed using a Finite Volume approach on a Cartesian grid. The 
expression of numerical fluxes does not affect the general coupling algorithm 
and we use a one-step high-order scheme proposed by Daru and Tenaud [Daru 
V,Tenaud C, J. Comput. Phys. 2004]. The Embedded Boundary method is 
used to integrate the presence of a solid boundary in the fluid. The coupling 
algorithm is totally explicit and ensures exact mass conservation and a balance 
of momentum and energy between the fluid and the solid. It is shown that the 
scheme preserves uniform movement of both fluid and solid and introduces no 
numerical boundary roughness. The efSciency of the method is demonstrated 
on challenging one- and two-dimensional benchmarks. 

Keywords: Compressible flows; Shock capturing scheme; Discrete Element 
method; Fluid-structure interaction; Embedded Boundary method 



1. Introduction 

This work is devoted to the development of a coupling method for fluid- 
structure interaction in the compressible case. We intend to simulate transient 
dynamics problems, such as the impact of shock waves onto a structure, with 
possible fracturing causing the ultimate breaking of the structure. An inviscid 
fluid flow model is considered, being convenient for treating such short time 
scale phenomena. The simulation of fluid-structure interaction problems is often 
computationally challenging due to the generally different numerical methods 
used for solids and fluids and the instability that may occur when coupling 
these methods. Monolithic methods have been employed, using an Eulerian 
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formulation for both the soUd and the fluid (for instance, the diffusive interface 
method [Hill]), or a Lagrangian formulation for both the fluid and the solid 
(for example, the PFEM method [26]), but in general, most solid solvers use 
Lagrangian formulations and fluid solvers use Eulerian formulations. In this 
paper we consider the coupling of a Lagrangian solid solver with an Eulerian 
fluid solver. 

A possible choice is to deform the fluid domain in order to follow the move- 
ment of the solid boundary: the Arbitrary Lagrangian-Eulerian (ALE) method 
has been developed to that end. It has been widely used for incompressible 
[TTl [TO] and compressible [K] fluid-structure interaction. However, when solid 
impact or fracture occur, ALE methods are faced with a change of topology in 
the fluid domain that requires remeshing and projection of the fluid state on the 
new mesh, which are costly and error prone procedures. Moreover remeshing is 
poorly adapted to load balancing for parallel computations. 

In order to allow for easier fracturing of the solid, we instead choose a method 
based on fictitious domains that solves the fluid flow on a fixed Eulerian mesh, 
on which a Lagrangian solid body is superimposed. A special treatment is then 
applied on fluid cells near the boundary and inside the solid. Different types of 
fictitious domain methods have been developed over the last thirty years. They 
can roughly be classified in three main classes: penalization methods, interpo- 
lation methods and conservative methods. Among penalization methods, the 
Immersed Boundary method is certainly the best known and most widely used 
for fluid-structure interaction. It was originally introduced by Peskin for incom- 
pressible blood flows [m [45] . The solid boundaries deform under the action of 
the fluid velocity, and the presence of the solid adds forces in the fluid formula- 
tion that enforce the impermeability of the solid. However, Xu and Wang have 
pointed out some numerical leaking of fluid into the solid [IH] . Following Lev- 
eque and Li [321 131] , they advocate the use of the Immersed Interface method, 
which incorporates jump conditions in the finite differences used. However, the 
absence of fluid mass loss is still not ensured exactly. In a different approach, 
Olovsson et al. [lOl |5] couple an Eulerian and a Lagrangian method by penal- 
izing the penetration of the solid into the fluid by a damped spring force. As 
the stiffness of the spring goes to infinity, the penetration goes to zero. Boiron 
et al. [S] and Paccou et al. [H] consider the solid as a porous medium, using 
a Brinkman porosity model. As the porosity goes to zero, the solid becomes 
impermeable. However, in both cases, as the stiffness grows or the porosity 
decreases, the use of implicit schemes is mandatory to avoid the severe stability 
condition of explicit schemes [S] I41j . For the high speed phenomena we con- 
sider, we use explicit solid and fluid solvers and an explicit coupling algorithm 
is better suited in order to avoid costly iterative procedures. 

Another class of fictitious domain methods consists in enforcing the bound- 
ary conditions through interpolations in the vicinity of the boundary, using the 
exact values taken by the fluid on the boundary [371 [13] . The method seems to be 
very versatile, being used with incompressible Navier-Stokes [371 [TB], Reynolds- 
averaged Navier-Stokes [lOl |29l [30] , turbulent boundary layer laws [6] and com- 
pressible Navier-Stokes [42]. The Ghost Fluid method developed by Fedkiw et 
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al. [m [T7] relies on the same type of principle for compressible fluids. The in- 
terface is tracked using a level-set function, and conditions are applied on both 
sides of the interface to interpolate the boundary conditions. The advantage of 
these methods is that they do not suffer from additional time-step restriction 
due to stability, and the order of accuracy of the boundary conditions can be set 
a priori. However, the interpolation does not ensure the conservation of mass, 
momentum and energy in the system. This can cause problems when dealing 
with shock waves interacting with solids. 

Conservative fictitious domain methods, generally referred to as Embedded 
Boundary methods, rely on a modified integration of the numerical fluid fluxes 
in the cells cut by the solid boundary. The original idea of the method can 
be traced back to Noh's CEL code [S^]. However, the approach brings out the 
problem of small cut-cells, where stability would require reducing the time step 
to very small values. Several solutions have been proposed to avoid time step 
restrictions in those cells, involving a conservative averaging in the vicinity of 
the boundary 03 [H Hi] . 

In this article, we develop a new coupling method for fluid-structure interac- 
tion within the framework of conservative fictitious domain methods. The finite 
volume scheme for the fluid calculation is modified in the near interface re- 
gion. The method can be implemented independently from the time integration 
scheme used for the fluid, whether based on space-time splitting or multi-level 
time integration. Conservative flctitious domain methods have proven to give 
satisfactory conservation results for inviscid compressible flows in the case of 
static solid boundaries. Nevertheless, to our knowledge, conservation issues of 
the coupling have not been studied in the case of moving solids. We establish 
new conservation results in such a case. Our coupling method is designed to be 
capable of treating the general case of moving deformable bodies. In the present 
work, however, we only consider non-deformable (rigid) solid bodies. The case 
of deformable bodies is the object of ongoing work. 

The fluid and solid solvers that we consider were chosen according to their 
ability to deal with shock waves and fracturing solids. The solid solver is based 
on a Discrete Element method, implemented in a code named MkaSD in the 
CEA [3S]. It can handle elasticity as well as fracture and impact of solids. Solids 
are discretized into polyhedral particles, which interact through well-designed 
forces and torques. The particles have a rigid-body motion, and fracture is 
treated in a straightforward way by removing the physical cohesion between 
particles. The work reported in this article is a first step towards the coupling 
with the MkaSD code. The time integration scheme used by Mka3D (Verlet for 
displacement of the center of mass and RATTLE for rotation (SHj) is retained for 
the rigid body treatment. Concerning the fluid solver, we use a Cartesian grid 
explicit flnite volume method, based on the high-order one-step monotonicity- 
preserving scheme developed in [5] and space time splitting. However we empha- 
size that our coupling method is independent from both the Discrete Element 
method (as long as a solid interface is defined) and the numerical scheme used 
for the fluid calculation. 

The article is organized as follows: we first present briefly the solid and fluid 
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methods in section [2] In sections [3] and [4] we describe the proposed exphcit 
couphng procedure between the fluid and the moving sohd in the framework of 
an Embedded Boundary method. The analysis of the conservation properties 
of the coupling is reported in section [5j where we show that mass, momentum 
and energy of the solid-fluid system are exactly preserved. In section |6j we 
demonstrate results about the preservation on a discrete level of two solid-fluid 
systems in uniform movement. Finally, we illustrate the efficiency and accuracy 
of the method on one and two-dimensional static and dynamic benchmarks in 
section [71 



2. Solid and fluid discretization methods 

2.1. Solid time-discretization method 

We consider a non-deformable solid (rigid body). The position and velocity 
of the solid are given, respectively, by the position of its center of mass X, 
the rotation matrix Q, the velocity of the center of mass V and the angular 
momentum matrix P. The physical characteristics of the solid are its mass m 
and its matrix of inertia R which, in the inertial frame, is a diagonal matrix 
with the principal moments of inertia /i, I2 and I3 on the diagonal. Here, we 
instead use the diagonal matrix D = diag((ii, ^2, da), where: 

Vze{l,2,3},d, = ^^-tif±i^-/,. 



The angular momentum matrix P can be related to the usual angular velocity 
vector J7 by the relation P = Dj(r2)Q, where the map j : E'' M'^^'^ is defined 
such that: 

Vx e M^ Vy e M^ j(x) • y = x A y. 

Let us denote by F and Ad the external forces and torques acting on the 
solid, and by At the time-step. In order to preserve the energy of the solid over 
time integration of the movement, we choose a symplectic second-order scheme 
for constrained Hamiltonian systems, the RATTLE scheme [53]: 

V"+5 = V" + — F", (1) 
X"+i = X" + AtV"+5, (2) 

pn+i ^ pn ^ ^j(_yvt»)Q« + ^A"Q", (3) 
Qn+1 ^ Qn ^ AiP"+5D-\ (4) 

with A" such that (Q"+i)'^Q"+i = 1, (5) 
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with A"+i such that (Q"+i)'^P"+iD-i + D-i(P"+i)'^Q"+^ = 0. (8) 

The symmetric matrices A" and A"+^ play the role of Lagrange multipliers for 
the constraints on matrices Q""*"^ and P""*"^. 

The scheme makes use of the velocity at half time-step V"+2^ which is 
constant during the time-step. Let us now consider the angular velocity. For a 
rigid solid, we have for all points x: 

X - X = Q ■ (XO - xO) , 

X'' and x*^ being material points of the solid at initial time. Using the identity 
ft A (Qx) = PD^^x for all x, the velocity at point x can be written as: 

V(x) = V + PD-i • (X" - x°) 

which is more convenient for use in the time scheme. In analogy with displace- 
ment, we consider P"+2 as constant during the time-step, and we define the 
velocity of point x at half time-step {n + 

V"+5(x) = V"+5 + P"+3D-i • (X" - x°) . 

2.2. Fluid discretization method 

The problem of the interaction of shock waves with solid surfaces can be at 
first studied using an inviscid fluid model. In this work, we consider inviscid 
compressible flows, which follow the Euler equations: 

wt -I- V • f(w) = 0, 

where w = {p, pu, pE)^ is the vector of the conservative variables, and f(w) is 
the Euler flux: 

/ pu 

f = pu(^u + p\ 
\ {pE+p)u 

where the pressure p is given by a perfect gas law: p = (7 — 1) (^pE — ^pu • u). 

To solve these equations, we use the OSMP numerical scheme, which is a 
one-step high-order scheme developed in |H1 H] . It is derived using a coupled 
space-time Lax-Wendroff approach, where the formal order of accuracy in the 
scalar case can be set at arbitrary order (in this paper, we use order 11, that is 
the OSMP 11 scheme). Imposing the MP conditions (Monotonicity Preserving) 
prevents the appearance of numerical oscillations in the vicinity of discontinu- 
ities while simultaneously avoiding the numerical diffusion of extrema. In one 
space dimension, on a uniform mesh with step-size Ax, at order p, it can be 
written: 
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where f^, , is the nth-order- accurate numerical flux of the scheme at the ceh 

interface (i + 5)- Given the I eigenvectors of the Jacobian matrix of the flux, 
the general expression of the numerical fluxes can be written: 

k 

where f is the standard Roe flux and the are corrective terms to obtain 
order p. Their exact expression can be found in [5]. At order the stencil of 
the scheme uses p + 2 points. 

In two dimensions, the fluxes are computed using a directional Strang split- 
ting which is second-order accurate. However the error of the scheme remains 
very low. 

3. Treatment of the cells cut by the solid boundary 

In order to take into account the position of the solid in the fluid domain, we 
rely on the Embedded Boundary method, which consists in modifying the fluid 
fluxes in cells that are cut by the solid boundary (named cut cells) , as in [23 [13] . 
At time i, for a cut cell C, we assume that the solid occupies a volume fraction 
ac- We also assume that the density, velocity and pressure are constant in the 
cell. The fluid mass, momentum and energy quantities contained in the cell are 
therefore equal to their value at the center of the cell times the volume of the cell 
and the volume fraction of fluid 1 — ac- In the same way, the computed fluxes 
are assumed to be constant on the faces of a cell. Denoting by KC1C2 the solid 
surface fraction of the face between cells Ci and C2 , the effective flux between 
Ci and C2 is the computed flux times the surface of their interface times the 
fluid surface fraction 1 — kc^^c^ ■ Additional fluxes come from the presence of the 
moving solid boundary. They are expressed in order to yield exact conservation 
of fluid mass and of the total momentum and energy of the system. 

For the sake of simplicity, we limit ourselves to two space dimensions. How- 
ever the three-dimensional case can be carried out in a similar way. Let us con- 
sider a fluid cell C cut by the boundary, as shown in Figure [Tj The indices /, r, t 
and b indicate respectively left, right, top and bottom in the sequel. Integrating 
the Euler equations on the cut cell and over the time interval [nAt, [n + l)At], 
and applying the divergence theorem, we get: 

,l-„-,A.,=A.(i^/„-i^/,. + i^/„ 

where Awq = the time increment and all fluxes are time-averaged 

over the time interval (the time averaging will be specifled later). At the solid 
walls, pressure forces cause momentum and energy exchange between the solid 
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Figure 1: Physical fluxes in a cut ceU 



and the fluid. They are taken into account through the exchange term Xjr. 
The detailed expression of Xjr will be given in section |T4j Finally, the quantity 
Awjr represents the amount of w" swept by each solid boundary present in 
the cell during the time step. Its expression will be given in section [4. 3[ 



4. Coupling algorithm 

Since the Discrete Element method is computationally expensive, the cou- 
pling algorithm should be explicit in order to avoid costly iterative procedures. 
In fact, the CFL condition of the explicit time-scheme gives the appropriate 
criterion for the capture of the high-frequency eigenmodes involved in the solid 
body fast dynamics. Moreover, as it is well known, explicit methods are more 
robust for impact problems. We choose the general structure of the algorithm 
as follows: 

— The position of the solid and the density, velocity and pressure of the fluid 
are known at time t 

— The fluid exerts a pressure force on the solid boundaries: knowing the 
total forces applied on the solid, the position of the solid is advanced to 
time t + At 

— The density, velocity and pressure of the fluid are then computed at time 
t + At. This step takes into account the new position and velocity of 
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the solid boundary, as well as the work of the forces of pressure on the 
boundary during the time step. 

The choice of the coupling algorithm is guided by the conservation of the 
global momentum and energy of the system and the conservation of constant 
flows (see section [6]). 

At the beginning of a time step, at time nAt, the position and rotation of 
the solid particle (X", Q"), the velocity and angular velocity of the solid particle 
(V",r2") and the fluid state w" are known. We choose the following general 
architecture for the algorithm: 



[ SOLID ] [ COUPLING ] [ FLUID ] 



X" Q'^ \'^'^ p" 



(3) Solid step 
(using predicted 
boundary pressure) 



(1) Computation 
of fluxes / 



(2) Predicted pressure is 
transferred to the solid 
boundary 



-| P.^^ Py. 



_J (4) Update of the boundary 
position, computation of the 
a"+i and k"+^ 



(5) Fluid update: 



p"+i = p" + Ap 
p"+iu"+i = p"u" + A(pu) 
p"+i_E"+i = |0"£;" + A(p_E) 



Steps (1) to (5) of the algorithm are computed successively, and are detailed 
in the following subsections. 

4-1. Computation of fluid fluxes and of the boundary pressure (steps (1) and 
(2)) 

Step (1) is a precomputation of fluxes without considering the presence of a 
solid boundary. As said above, the fluxes are computed in every cell using the 
OSMPll scheme. However, we emphasize that the coupling algorithm does not 
depend on the choice of the numerical scheme. The fluxes are then stored for 
later use in step (5). 

The other aim of this step is the computation of mean pressures in each 
cut-cell during the time-step in each direction and Py. These pressures, 
transferred to the solid boundary in step (2), account for the forces exerted 
by the fluid on the solid during the time-step. The same mean pressures will 
be used in step (5) to compute the momentum and energy exchanged between 
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the solid and the fluid. In this way, the choice of and Py has no effect on 
the conservation of fluid mass, momentum or energy of the system. On the 
contrary it is a key ingredient for the exact conservation of constant flows (see 
section |6| . The explicit structure of our solid and fluid methods allows several 
possibilities for the choice of boundary pressures while maintaining the stability 
of the coupling algorithm. This is unusual in fluid-structure interaction. 

For the Strang splitting algorithm used here, the mean pressures p^ and 
Py are the pressures used for computing the fluxes in the x and y directions, 
respectively. An analogous relation could be derived for other time integration 
methods, such as Runge-Kutta for instance. 



4-2. Computation of the solid step (step (3)) 

Step (3) consists mainly in the application of the time integration scheme 
for the rigid body motion described in section [ZT] The essential difference with 
an uncoupled version lies in the integration of boundary pressure forces. As we 
consider an explicit coupling, the only boundary pressures available are p^ and 
Pv 




Figure 2: Geometric description of tlie particles 



The solid is assumed to be polygonal (in two space dimensions) as described 
in Figure [2] We denote by ^ the list of all faces of the solid in contact with 
fluid. For every face J- & ^, the position of the center of the face is given by 
vector X_F, and we denote by Sjr its surface and n_F its normal vector (oriented 
from the solid to the fluid). The fluid pressure force Fjr exerted on face £ ^ 
is written as: 

• e, = -p^S^n^^ (10) 

• = -PyS^n^^ (11) 

The total fluid pressure force is the sum of the contributions on each face: 

F'; - E (12) 
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The fluid pressure torque Ai^ is the sum of the torques of the pressure forces 
at the center of mass of the soUd body: 

M] = ^ A (X" - X^) 

The sohd time-step is written as in equations ([T]) to ([8|, with the only dif- 
ference that the fluid pressure force and torque are taken constant during the 
whole time-step, equal to and (including in equations (j6| and (j7|). The 
fact that FJ, V"+3 and P"-+5 are constant during the time-step will be 

used in the conservation analysis in section [Sj 

4-3. Update of the boundary and of the volume fractions (step (4)) 

Several tasks are carried out in step (4). For each cell C, the new solid 
volume fraction of the cell a^"*"^ surface fractions n^^^ are computed. 

In addition, for each solid boundary the pressures p,j. and Py are stored and 
the swept quantities Aw^ used in ^ are evaluated. 

In two dimensions, the solid boundary is polygonal, and we therefore only 
have to deal with plane boundaries J-. In order to simplify the computation of 
the average of and Py on we also assume that each boundary is contained 
only in one cell at time nAt. The computation of the contribution of Awjr to 
each cell is also easier if T is entirely in the cell at time (n -I- l)At. We denote 
4>„(J") the position of boundary J" at time nAt. We subdivide the polygonal 
boundary of the solid into a set of plane boundaries J" such that each of $„ (J") 
and i>„+i(.7-") is contained in one cell (which may differ, as shown in Figure [s]). 
Each newly computed boundary € ^ stores every variable necessary for the 
coupling: the surface Sjr and the normal vector njr of $„(J-"), the center of mass 
Xjr of T, and we define X^r = $o(Xjf')- The boundary also stores the pressures 
p^ and Py in the cell occupied by $„(J"), and the velocity of the center of the 

boundary at time (n + ^)At, Vjr ^ , computed as: 

Y»+i ^ yn+i _^ pn+iQ-l . - X^r) (13) 

The swept quantities Au;^ are computed as the integral of w" in the quad- 
rangle bounded by $,1(7^) and $„+i(J^) (see Fig. [s]). The condition 

J2^w'^-J2k^'-^c)w^ (14) 

c 

is then automatically satisfied as the set of such quadrangles is a partition of 
the volume swept by the solid during the time step. 

The computation of a^^"*"^ and involves the intersection of planes with 
rectangles, and can be carried out geometrically. The application of the diver- 
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Solid boundary at t = nAt 



Figure 3: Update of the boundary and computation of the Aw'^ 

gence theorem to ccU C shows that the following relations are satisfied: 



(15) 
(16) 



These conditions correspond to the geometric conservation laws (GCL) in ALE 
methods [32 , and will be used in the analysis of the consistency of our method. 

4-. 4- Modification of fluxes (step (5)) 

Step (5) of the algorithm consists mainly in the computation of the final 
values w^'^^ in each cell, using a fully discrete expression of Eq. It is the 
only part of the algorithm where the fiuid "sees" the presence of a solid. The 
explicit fluid fluxes were pre-computed in step (1) on the Cartesian regular grid, 
and the modification of fluxes aims at conserving the mass of fluid and balancing 
the momentum and energy transferred to the solid during the time-step. 

The exchange term in ^ can be written as 

Tec 

where fjr is the fluid flux at the solid boundary (see fig. [T]) , that are approx- 
imated as: 



S 



(17) 



Here 'Vjr ^ is the velocity of the center of the boundary and is defined in ( 13 ) 
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Using the fluid fluxes given by the OSMPll scheme, we flnafly compute the 
time increment Awc from the foUowing fuhy discrete version of equation ^ : 



The value of wc is then updated in every cefl: wji'^^ = Wq + Awc- 

A main difference with ^14; hes in the time integration of ceU-face apertures 
(1 — Kc)- Falcovitz et al. [13] use time-averaged cefl-face apertures over the time 
step (at time (n + 5) At), ensuring consistency (in the sense that the uniform 
motion of a sohd-fluid system is exactly preserved). Here we instead take 
but the analysis of the coupling outlined below does not depend on the chosen 
time for k and would still hold if we took k"+2 instead. This theoretically re- 
duces the accuracy of the method in cut-cells, however our numerical results did 
not justify the added complexity of computing the time average k"+3. More- 
over, using k"'^^ , we still ensure consistency but also additional properties as 
shown in section [6l 

In order to avoid the classical restriction of the time-step due to vanishing 
volumes: 

(1 - ac)i'nm{Ax,Ay) 
i|ui|+c 

where c is the local speed of sound, we resort to the mixing of small cut cells 
with their neighbors to prevent instabilities. 

4.5. Conservative mixing of small cut cells 

Two main methods have been developed to ensure the stability of conser- 
vative Embedded Boundary methods. A first method consists in computing a 
reference state using nonconservative interpolations, modified by redistributing 
the conservation error on neighbouring cells |43l [121 IMj • A second method is to 



compute a fully conservative state using a formula similar to (18). For stabil- 
ity reasons, small cells are merged with neighboring cells using a conservative 
procedure (originating from Glimm's idea [H]). We choose this second class of 
method, as ^25^ and (14j . 

Target cells need to be defined for small cells to be merged with them. [25] 
defines an equivalent normal vector to the boundary in the cell and mixes the 
cells preferentially in that direction. jl4| rather merges newly exposed or newly 
covered cells with full neighbours having a face in common. In order to deal with 
cells occupied by several boundaries (impact of two solids), we cannot define a 
normal vector in every cell and we choose to improve the strategy applied in 
[H]. We define small cells as ac > 0.5. For mixing two cells C and Ct, so they 
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have equal final value w, the following quantities are exchanged: 



M, 



CtC 



ac + o^Ct 
ac + o-Ct 



{wct - wc) 
{wc ~ wcj.) 



and it is easy to check that wc + Mcct — '^Ct + ^'^CtC- In the two dimensional 
case, we select the target cell Ct as the fully-fluid cell {acj, — 0) nearest to 
C, such that the path between the two cells does not cross a solid boundary. 
A recursive subroutine finds such a target cell in a small number of iterations, 
without any restriction on the geometry of the fluid domain. 

5. Analysis of the conservation of mass, momentum and energy 

In this section, we analyze the conservation properties of the coupling algo- 
rithm. These properties are verified for periodic boundary conditions or for an 
infinite domain. 



5.1. Integration on the fluid domain 

Integrating w on the fluid domain at time {n + l)At, we obtain using 



( 18 ) and the cancellation of fiuxes on each cell face: 



A 



:5:(l-ar>c+E^/^ 



c 



Using (14) we finally get: 
1 



AxAy Jf2^+i 



w 



ji+i 



AxAy' 



AtS 



E(i-cK + Ea^/- 



AxAy 



At 



AxAy ^ 



(19) 



the expression of /jr being given in Eq. [TT] 



The first component of system (19 1 expresses the fluid mass conservation. 



In order to proceed with the analysis of momentum and energy conservation, 
let us now turn to the solid part. 

5.2. Solid conservation balance 

Since the solid is treated using a Lagrangian method, the conservation of 
solid mass is straightforward. The fluid pressure force applied on the solid during 



the time step is given by (10 1, (11) and (12). Let us consider a solid boundary 



G 5", and denote by A'Pjr the solid momentum variation induced by the 
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pressure forces on T , and jf the corresponding energy variation. Recalling 
that the pressure forces are kept constant during the time step, the balance of 
momentum and energy is given by: 

Af^ = A^F^ ■ (i^ X ^^^'^^ " ^^^^ ' ^'^^'^ 



Finally, using the expression of forces Fjr, we obtain: 
AP^ = -At^p,n^ 
AVI = / Vy^^T 



Comparing with section |5.1[ the balance of momentum and energy in the 
fluid domain results in: 



j p"+i£;"+i + V A5^ = ( p'^E" 



This demonstrates the conservation of momentum and energy for the coupled 
system. 



6. Conservation of constant flows 

In this section we analyze the consistency of the coupling method, in the 
sense defined in }14j . meaning exact conservation of uniform flows by the cou- 
pling algorithm. Two cases are analyzed. The flrst one, also considered in |14j . 
consists of a solid immersed in a fluid and moving at the same velocity. The 
second one, not considered before, demonstrates the correct representation of 
the slip boundary condition along walls. These simple cases have been a guide 
to design the algorithm, as the preservation of such flows is a basic criterion for 
the quality of the method. 

In the whole section, we consider a constant fluid state: p" = pq, u" = uq, 
v" = vq and ~ po everywhere. The fluxes / are such that = fi — 
{poUo,poul + PcPoWqWo, (poco + Po)uo)'^ and ft = ft, = (poVq, poUqVo, pav^ + 
Po, (poco + Po)uo)'^ . In this case, the efficient pressures on the boundary of the 
solid are p^ —Py = Po- 
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6.1. Steady constant flow with moving boundaries 

We consider an arbitrarily shaped rigid body, moving at constant velocity 
with no rotation, immersed in a uniform fluid flowing at the same velocity. 

The solid is a closed set, and we denote by fi" the solid domain at initial 
time. We have: 



V 5^n^ = (f ndS^O 



Using (10) and (11), we obtain: 

This induces: 
In the same way, 

V = -J^poS^n^ A (Xr - X^) ^po(f (X^ - X) A nd5 = 
J, J, Jan^ 

The volume swept by boundary J- is AtSjr{uo, vq)^ ■ njr. Since the initial state 
is constant, Awjr is given by: 

ATT ^ ^ / tj \ 

AxAy 

In addition, as the solid translates without rotation, the normal vector njr to 



a boundary J- is constant in time. Using this property in equations ( 15 ) and 



( 16 ), we easily conclude that (1 — ac)Awc = 0. Thus = w", showing that 

the constant flow is left unchanged by the coupling algorithm and the mixing 
of small cells. 

6.2. Free slip along a straight boundary 

We consider an undeformable, fixed solid consisting in a semi-infinite half- 
space. The solid boundary is a straight planar boundary with a constant normal 
vector n such that: 

uo • n = (20) 

This initial state describes the free slip of the fiuid along the straight boundary. 
In the inviscid case, no boundary layer should develop in the vicinity of the 
boundary. The conservation of such fiows ensures that the boundary is not seen 
by the fluid as being artiflcially rough. 

As the solid is fixed, ac and kc remain constant over time and Aw^ is equal 



to zero. From equation (18), and using (15), (16) and (20), the components of 
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Awe are calculated as: 



(1 - ac)Apc = -At 

(1 - ac)A{pv)c = -At J2 
Tec 

{l-ac)A{pE)c^-AtY, 
Tec 

This shows that the constant flow is preserved by step (5) of the algorithm. 
This result is not modified by the mixing procedure. We thus have shown the 
exact preservation of the free slip of the fluid along a straight boundary. 

7. Numerical examples 

In the following, we consider a perfect gas, with 7 = 1.4. In all computations 
the CFL number was fixed equal to 0.5. 

1.1. One- dimensional results 

A piston of density 2 kg.m~'^ and length 0.5 m is initially centered at x = 
2 m, in a one-dimensional, 7m-long tube, whose ends are connected by periodic 
boundary conditions which allow an easier comparison with ALE results. The 
gas initial pressure and density are equal to 10^ Pa and 10 kg.m~^ for x < 2m 
and X > 5m and to 10^ Pa and 1 kg.m"'^ elsewhere. The system is initially at 
rest. The initial pressure difference between the two sides of the piston triggers 
its movement and the propagation of waves in the fluid regions (a rarefaction 
in the left region and a shock wave in the right one). Wave interactions then 
occur at later time. The fluid pressure at time t = 0.003s is shown in Figure |4] 
and the trajectory of the solid is presented in Figure [5j The x — t diagram over 
a longer time (0.01 s) is shown in Figure [6j 

An ALE computation was done for comparison, using a uniform grid mov- 
ing at the solid velocity. The solid position and velocity are updated using 
the same second-order Verlet scheme. We compared the numerical results ob- 
tained through the Embedded Boundary method on 100, 200, 400, 800, 1600, 
3200, 6400 and 12800 points grids with a 51200-points ALE grid, considered 
as the reference solution. We observe a second-order convergence of the solid 
position (Figure [7| and a super-linear convergence of order 1.2 of the fluid pres- 
sure (Figure |8]). The convergence rate is optimal for the solid (Verlet scheme 
is second-order accurate). The convergence rate for the fluid pressure is not 
optimal, due to the presence of discontinuities, but is not affected by the solid 
coupling. 



St 
AxAy 

S 



n • Ug = 



4- ((n ■ Uo)mo + Pon^) + ^ ^^^^ Pon^ = 



AxAy 



Tec 



AxAy 



s 



AxAy 



4-(n- Uo)(poeo +Po) = 
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Figure 4: Pressure at time t = 0.003 s 
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Figure 5: Time evolution of the solid position 

7.2. Double Mach reflection 

A Mach 10 planar shock wave reflects on a fixed 30° wedge, creating a Mach 
front, a reflected shock wave, and a contact discontinuity which develops into a 
jet along the solid boundary. This benchmark was first simulated on a Cartesian 
grid aligned with the solid boundary, using different finite volume methods 
[27| l47l [8] . Non-aligned grid methods were also tested on this benchmark, using 
Embedded Boundary methods |43|, [7], non-conservative Immersed Boundary 
methods [5T], h-hox methods [53], and kinetic schemes 31i. The position of 
the tip of the jet is an important characteristic of the accuracy of the results. 
The comparison with grid-aligned results shows that it is better recovered by 
conservative methods than non- conservative methods [43l [7] . 

We have simulated the problem on a grid aligned with the wedge (aligned 
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Figure 6: x — t diagram (the position of the solid is in deep blue) 
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Figure 7: Convergence of the solid position L°°-error 



case, Figure 9| and on a grid aligned with the incident shock wave (non-aligned 
case, Figure 10). The two results are very similar, and agree with [7, 43, 27] |H]- 



One can remark that all the features of the flow are captured at the correct 
position in the non-aligned case. The jet propagates along the wall without 
numerical friction due to the conservation of free slip along a straight bound- 



ary (section 6.2). In the principal Mach stem, the discontinuities are slightly 
more oscillatory than in the aligned case. This can be identified as a post-shock 
oscillation phenomenon to which Roe's scheme is especially prone (see, for in- 
stance, [2H1I1]), and is not related to the coupling method. Nevertheless, the 
perturbations stay localized in the vicinity of discontinuities. 
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Figure 8: Convergence of the fluid pressure L^-error 




Figure 9: Aligned case: 30 contours of fluid density from 1.73 to 21, Ax = Ay = 1/220 



7.3. Lift-off of a cylinder 

This moving body test case was first proposed in fT^, using a conservative 
method. A rigid cylinder of density 7.6kg.m^'^ and diameter 0.1 m, initially 
resting on the lower wall of a Im x 0.2m two-dimensional channel filled with air 
at standard conditions (p = 1 kg.m~"^, p — 1 Pa), is driven and lifted upwards 
by a Mach 3 shock wave. Gravity is not taken into account. The problem was 
simulated in [201 ISl US]- In Figure 11 we present our results on a uniform 
1600 X 320 grid at times 0.14 s and 0.255 s. The cylinder was approximated by 
a polygon with 1240 faces. 

Our results agree well on the position of the solid and of the shocks with 
Arienti et al. [7 and Hu et al. [25] . However, some differences should be noted. 
First of all, some reflected shock waves in our results seem to lag slightly behind 
their position in previous results. This difference might be caused by small dif- 
ferences in the final position of the solid. Hu et al. |2S] also discuss the presence 
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Figure 10: Non-aligned grid case: 30 contours of fluid density from 1.73 to 21, Ax = Ay = 
1/220 

of a strong vortex under the cylinder in the results of Forrer and Berger j20| . 
They dismiss it as an effect of the space-time splitting scheme employed which 
affects the numerical dissipation. We also obtain this vortex, which does not 
disappear as we refine the mesh. We rather believe that this vortex is associated 
with a Kelvin-Helmholtz instability of the contact discontinuity present under 
the cylinder (Figure [T2|). 

In Figures [13] and |14| we present convergence results on the final position 
of the center of mass of the cylinder, compared to those of Hu et al. [25] ■ We 
observe that our results exhibit a fast convergence process, which is not the case 
in [251 . Let us note that no exact solution exists for the final position of the 
cylinder. The final position we found is however in the same range as in |25j . 
The results also compare well with Arienti et al. [3]. The improvement lies in 
the combination of the conservative interface method |25j with a conservative 
coupling and a second-order time-scheme for the rigid body motion. For this 
difficult case, the maximal conservation relative errors due to coupling were 
bounded by 4.10"^ over the whole simulation time, and no drift was observed. 

In Figure [15] we present the relative computational cost of the coupling. The 
relative cost is defined as the ratio of the computational times dedicated to the 
coupling method and to the fiuid and solid methods. In the rigid body case, the 
cost of the solid method is very low compared to that of the fluid method. As 
the coupling method is explicit and local, the computational cost is located on 
a manifold one dimension lower than the dimension of the whole space. In the 
two-dimensional case, the coupling is on a one-dimensional manifold. Indeed, we 
observe in Figure [15] that the relative cost of the coupling decreases as the grid 
is reflned, with a slope of 0.5, and that the coupling cost remains lower than the 
fluid and solid costs, amounting to approximately 10-20% for the grids yielding 
sufficient accuracy. 
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Figure 12: 60 contours of fluid density from to 12 at final time, Ax = Ay = 6.25 X 10 ^ 



7.4- Flapping doors 

We propose this new fluid-structure interaction case as a demonstration of 
the robustness of our approach and as a first step towards fracture and impact 
simulations. The flapping doors case involves separating or closing solid bound- 
aries, with cells including several moving boundaries. The algorithm is shown to 
be able to deal with such difficulties. Two doors initially close a canal and are 
impacted from the left by a Mach 3 shock. The canal consists of two fixed rigid 
walls, 2m long and 0.5m apart. Each door consists of a 0.2-m long and 0.05m- 
wide rectangle, completed at both ends by a half-circle of diameter 0.05m. The 
doors are respectively fixed on points (0.5,0.025) and (0.5,0.475) at the center 
of the half-circles. They can rotate freely around these points. The Mach 3 
shock is initially located at a; = 0.43m. The density of the solid is 0.1 kg.m^'^ 
and the pre- and post-shock state of the fluid are {p,u,v,p) = (1,0,0,1) and 
{p,u,v,p) = (3.857,2.6929,0,10.333). In Figure [l6| we show the density field 
obtained using a 1600x400 grid, at times 0.125s, 0.25s, 0.375s and 0.5s. After 
the incident shock hits the doors, it reflects to the left and the doors open due 
to the high rise in pressure. The opening of the doors produces a jet preceded 
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Figure 13: Comparison of the convergence of the horizontal position of the center of mass of 
the cyhnder with Hu et ah I25| 
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Figure 14: Comparison of the convergence of the vertical position of the center of mass of the 
cylinder with Hu et al. j25j 

by a shock wave propagating to the right. Then complex interactions of waves 
occur, due to door movements and interaction with the waUs. Kelvin-Helmholtz 
instabihties of contact hues can be observed at t=0.5s. It is worth noting that 
symmetry of the flow about the centerline of the canal axis is remarkably well 
preserved by the coupling method. 

As the doors remain tangent to the canal walls during their rotation, the 
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Figure 15: Ratio between the coupling computation cost and the fluid and solid costs 

fluid cannot pass between the wall and a door at its hinge. When the doors 
approach the walls at maximum rotation, the fluid is compressed, and eventually 
pushes them back. This is observed at time t — 0.2162s and t — 0.486s in Figure 
[TTj which presents the time evolution of the doors rotation angle. In the first 
case, the distance between each door's straight boundary and the wall is less 
than 0.002m, while the size of a fluid cell is Ax = 0.00125m. The method is 
able to deal with the fact that most cells along the wall are cut by the moving 
boundary and contain several moving boundaries. Treating this test case with 
an ALE method would require several remeshings in the course of the simulation, 
especially in the initial separation of the tangent door tips, and when the doors 
approach the walls. 



8. Conclusion 

We have presented a new Embedded Boundary algorithm for coupling a Fi- 
nite Volume method for compressible fluid flows with a rigid body. This explicit 
algorithm has the advantage of preserving the usual CFL stability condition: the 
time-step can be taken as the minimum of the full cell size fluid and solid time- 
steps. The combination of the Embedded Boundary method for the flctitious 
fluid domain and of the coupling strategy ensures the conservation of fluid mass 
and the balance of momentum and energy between fluid and solid. In addition, 
the exact conservation of the two constant states described in section [6] gives 
good insight on the consistency of the method: we prove a form of Galilean in- 
variance and the fact that the treatment at the boundary introduces no spurious 
roughness or boundary layers. The numerical examples suggest the second-order 
convergence of the solid position and the super-linear convergence of the fluid 
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Figure 16: Density contours at times t = 0.125s (a), t = 0.25s (b), t = 0.375s (c) and t = 0.5s 
(d) 

state in norm, while our results on two-dimensional benchmarks agree very 
well with body-fitted methods and improve Immersed Boundary results. We are 
also capable of dealing with solid boundaries moving close to each other, which 
is promising for impact simulations. The method is computationally efficient, 
as the coupling adds an integration on a space one dimension smaller than the 
fluid and solid computation spaces. The present method is therefore perfectly 
liable to be extended to a deformable solid, and was designed to extend natu- 
rally in three space dimensions. The remaining difficulty is the ability to define 
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Figure 17: Time evolution of the rotation of a door 



and track the solid boundary surrounding a Discrete Element assembly. 
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